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Molecular Dynamics studies of chemical processes in solution are of great value in a wide spectrum of ap¬ 
plications, which range from nano-technology to pharmaceutical chemistry. However, these calculations are 
affected by severe finite-size effects, such as the solution being depleted as the chemical process proceeds, 
which influence the outcome of the simulations. To overcome these limitations, one must allow the system 
to exchange molecules with a macroscopic reservoir, thus sampling a Grand-Canonical ensemble. Despite the 
fact that different remedies have been proposed, this still represents a key challenge in molecular simulations. 

In the present work we propose the Constant Chemical Potential Molecular Dynamics (C^MD) method, 
which introduces an external force that controls the environment of the chemical process of interest. This 
external force, drawing molecules from a finite reservoir, maintains the chemical potential constant in the 
region where the process takes place. We have applied the C^MD method to the paradigmatic case of urea 
crystallization in aqueous solution. As a result, we have been able to study crystal growth dynamics under 
constant supersaturation conditions, and to extract growth rates and free-energy barriers. 

J. Chem. Phys. 142, IffllS (2015); http://dx.doi.org/10.1063/1.4917200 


I. INTRODUCTION 

Physical-chemical processes taking place in liquid solu¬ 
tion are ubiquitous, and are at the heart of a wide variety 
of applications. Electro-chemical reactions, surfactants 
adsorption, self-assembly, crystal nucleation and growth, 
are just a handful of such applications. However, in many 
cases, a detailed description of the molecular mechanisms 
at play during these processes is still lacking. 

Molecular Dynamics (MD) represents a powerful 
method for investigating such processes with atomistic 
detail. However MD suffers from many limitations such 
as the relatively small time and size scales that can be 
simulated. With the presently available computational 
resources, classical MD calculations, based on empiri¬ 
cal potentials, can typically study systems of size up to 
10'^ 10® atoms. Such size limitations are particularly 

dramatic in the simulation of phase transformations, as 
in the paradigmatic case of crystal growth from solution. 
In such a case, while the crystallization proceeds, the so¬ 
lution is depleted, thus changing its chemical potential 
and affecting the growth process itself. For this reason 
MD results require seizable finite-size corrections before 
being compared with the experimental results, which in¬ 
volve much larger system sizeJ^^t^. 

In general, the numerical approach to prevent such fi¬ 
nite size problems is that of sampling configurations in 
the Gran-Canonical (GC) ensemble, namely simulating 
a system in contact with an external reservoir, which 
maintains the chemical potential constant. However, the 
classical methods for GC samplin^^HH require on-the-fly 
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insertion and removal of particles. This constitutes a cru¬ 
cial obstacle for chemical processes in solution, since the 
probability of a successful particle insertion in a dense 
fluid is extremely lov^l, preventing efficient numerical 
studies. 

Recently, the development of adaptive resolution sim¬ 
ulation methods AdResS (s ee e.g . Ref. [ini and refer¬ 
ences therein) and iif-AdRes^ ^^b^ l ^ri alter¬ 

native path to MD simulation in the GC ensemble. In 
these methods the region of primary interest is described 
with atomistic detail, while its surroundings are treated 
through a lower resolution model. It has been shown 
that the low resolution region can act as an external 
molecule reservoir, enforcing the sampling of the GC en¬ 
semble in the atomistically resolved region. Moreover, in 
this low resolution region, particle-insertion or swapping 
techniques can be applied more efhcientljff^. 

Here, inspired by this approach, we use an analogous 
volume subdivision and implement a method that ad¬ 
dresses the problem of solution depletion in MD simula¬ 
tions. To be more definite, we focus our attention on 
the simulation of crystal growth from solution, which 
lately has received much attentioill^lHSI and for which 
the pitfalls of limited size numerical modeling have been 
underlinecP^. Nonetheless we stress that the present 
method can be applied to many other systems. 

In our scheme the region containing the growing crys¬ 
tal and its immediate surroundings is maintained at con¬ 
stant solution concentration, while the remainder of the 
simulation box acts as a molecular reservoir. In this way 
we are able to study the crystal as it grows or dissolves 
in contact with a solution at constant chemical potential. 

The paper is organized as follows: first, in Sec.|Hj our 
method is introduced and described in detail. After that 
we apply it to the crystallization of urea in aqueous so- 
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FIG. 1. Qualitative behavior of the solute local density in 
the vicinity of a planar crystal surface. The blue, vertical line 
indicates the solid-liquid interface position at z = zi. The 
light blue shaded area corresponds to the TR, dehned in the 
text. Within this region, the dashed line is explanatory and 
non-representative of the realistic Uu behavior. 


FIG. 2. Same proBle displayed in Fig. with the addiction 
of an external force applied at z = zj + Dj? (green vertical 
line). The green shaded area corresponds to the transition 
region connecting the GR and the reservoir. Also in this case 
the dashed profiles are purely explanatory. 


lution. The technical details of the calculations are pre¬ 
sented in Sec. |III| Then, in Sec.|IVl the simulation results 
are illustrated and discussed. Finally, Sec. |V]is devoted 
to the conclusions and possible perspectives of our work. 


II. CONSTANT CHEMICAL POTENTIAL MD 

In the present section we describe in detail our 
new method, named Constant Chemical Potential MD 
(C/iMD), which allows to perform MD of chemical pro¬ 
cesses in a solution at constant chemical potential. In 
order to illustrate the method we focus on the problem 
of crystal growth in a binary solution, while in princi¬ 
ple any first order phase transition can be considered. 
Moreover, for the sake of simplicity, we consider a planar 
crystal-solution interface. Extension to other geometries 
is possible, although rather more complicated. 

In Fig. we show a qualitative description of the lo¬ 
cal solute density Uu in the neighborhood of the crystal- 
solution interface, located at zj. To the left of the inter¬ 
face there is the crystal region, characterized by the solid 
density n^. On the other side of zi we find a Transition 
Region (TR), of length in which the growth process 
takes place. The extension of the TR is such that for 
z> zi+^ the density reaches its bulk value n®. The Uu 
profile within the TR depends on the kinetics of crystal¬ 
lization and on the diffusivity of the liquid, and cannot be 
easily guessed a priori. Here we are implicitly assuming 
the considered transition does not involve macroscopic 
correlation lengths, so that ^ is finite. 

During crystal growth or dissolution, the solid-liquid 
interface moves together with the TR. In a macroscopic 
system nP remains unchanged and the density at the 
boundaries of the TR is constant, leading to a stationary 
growth process. This is in contrast with the behavior of 
a finite-size system, in which nP varies due to the limited 
number of molecules. 

Our method aims at restoring this stationarity condi¬ 
tion, enforcing a constant solution composition in a Con¬ 
trol Region (CR) in contact with the TR, while the rest 
of the solution volume is used as a molecule reservoir. 


as shown in Fig. To control the solution density, an 
external force is applied at a fixed distance Dp from 
the moving crystal interface. acts as a membrane, 
regulating the flux of molecules between the CR and the 
reservoir, in order to maintain the former at a constant 
concentration. If we assume that the effect of the exter¬ 
nal force is felt in a region of size ap, the CR is located 
between zi -I- ^ and zi + Dp — apl2. In order not to 
affect the growth process the CR should be larger than 
the typical correlation lengths of the solution. We stress 
the fact that the force is applied at a fixed distance from 
the solid-liquid interface, to maintain a stationary solu¬ 
tion environment in the proximity of the growing crystal 
interface. 

As the crystallization proceeds the reservoir is de¬ 
pleted, so that the chemical potential can only be main¬ 
tained for a limited amount of time Tp. In general, when 
applying the C/xMD technique to some molecular pro¬ 
cess, this time limit must be taken into account, and the 
simulation box has to be properly engineered so that 
is able to act efficiently during the timescale of interest. 

We write the external force F^ as: 

Ftiz) = hinf^-no,)G{z,Zp), ( 1 ) 

where i labels the solute or solvent species, ki is a force 
constant, noi is a target density and G is a bell-shaped 
function which is nonzero in the vicinity of the force cen¬ 
ter Zp = zi + Dp. Eq. 0 defines a harmonic-like force, 
localized at a fixed distance Dp from the solid, which acts 
on the solution molecules to compensate the deviations 
of the instantaneous CR density nf^ from noi- 

If Ni is the total number of i-species molecules and 
the CR volume, we evaluate as: 


Ni 


n?^ = 


yCR 




1 = 1 


where: 


= 


if Zj G CR 
otherwise 


( 2 ) 


( 3 ) 
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is a function that selects the molecules inside the CR. In 
the practice, we let 9(z) switch continuously to 0 at the 
CR boundaries, in order to avoid sudden jumps in nf^. 

We then define the function G{z,Zf), localizing the 
action of close to the force center Zp, as: 


Gw{z 



1 + cosh 




( 4 ) 


is nonzero only for z ^ Zp, with an intensity peak 
proportional to w~^ and a width proportional to w. We 
observe that is non-conservative, since it is not pos¬ 
sible to define a potential function that leads to Eq. 0. 

As a further remark we note that Eq. 0 defines two 
separate forces acting on the solute and solvent species. 
In principle the concentration of both species affects 
the chemical potential of the solutioiP^ and should be 
controlled. However, we are mainly interested in con¬ 
stant pressure simulations, where only a single popula¬ 
tion needs to be subj ect to because the action of 
the barostat algorithrcPHSl] guarantees a prompt equili¬ 
bration of the other species concentration. If instead an 
NVT dynamics is considered, F^ should act on both so¬ 
lute and solvent species. 

To conclude this section let us summarize the C/rMD 
scheme in explicit algorithm steps: 

1. the solid-liquid interface position zi is located on- 
the-fly analyzing the solvent molecules distribution 
within the box, 

2. the force center Zi?, is updated maintaining a fixed 
distance Dp from zi, and the CR position is shifted 
accordingly, 

3. the densities are evaluated via Eq. 0. 

4. The MD equations of motions, including the exter¬ 
nal forces of Eq. 0 , are integrated. 

We recall once again that this method is effective for a 
finite validity time Tp, depending on the availability of 
molecules in the reservoir. 


III. TEST CASE AND SIMULATION SETUP 

In the present section we introduce the typical setup of 
our calculations. We considered here the case of a urea 
crystal growing in aqueous solution in slab geometry, as 
shown in Fig. The system is generated starting from 
an approximately 2.5 nm thick crystal slab, periodically 
repeated in the x and y directions. Along z, the slab 
exposes either the {001} or the {110} face, the two sta¬ 
ble faces in urea crystal growth from aqueous solutiorP^. 
Such a crystal is immersed in a supersaturated solution, 
generated by means of the PackmoP^and genbox (GRO- 
MACS package) utilities. The typical size of the simu¬ 
lation box is of 2.7 X 2.7 x 14nm^, Periodic Boundary 
Conditions (PBC) are applied in the x,y and z directions. 



FIG. 3. Typical simulation box for the study of urea crystal 
growth in aqueous solution. The vertical lines indicate the 
force center Zp, where F^ is applied (see Eq. 0), located 
at both sides of the slab, at distance Dp from the crystal 
interfaces. The gray-shaded areas indicate the CR. 


The number of urea and of water molecules for each sim¬ 
ulated system is reported in Tab. |lj The initial solute 
concentration has to be larger than the target concentra¬ 
tion, to allow a stationary growth regime for a timescale 
of 10 - 100 ns. 

Following our previous worlPsE^j_TO use the Gener¬ 
alized Amber Force Field (GAFFjpSPH to model urea, 
and the T1P3P modeP^ for water. The molecular bonds 
are constrained through the LINGS algorithnP^ and long 
range electrostatics is handled by means of Particle Mesh 
Ewald (PME) methocP^. 

The system is kept at a constant temperature of 300 
K and pressure of 1 bar by using the stochastic veloc¬ 
ity rescaling thermostalP^ and the Parrinello-Rahman 
barostat!^. Because of the planar symmetry, the barostat 
is applied in its semi-isotropic version, so that the z-side 
of the box, parallel to the growth direction, is rescaled 
separately from the x and y sides. Since the system is 
rescaled during the dynamics, also the G/rMD length pa¬ 
rameters, as e.g. Dp, have to be rescaled accordingly. 
Our implementation redefines the lengths on-the-fly, fol¬ 
lowing the fluctuations of the z-edge of the box. In the 
considered simulations F^ does not affect the pressure of 
the system by a significant amount, since it determines 
a surface effect, small compared to the potential and ki¬ 
netic bulk contributions. For this reason in the presented 
results we did not couple the external force with the baro¬ 
stat equations. 

After the initial equilibration phase, a 2 fs integration 
step has been chosen for all the production runs. All 
the calculations have been performed using the GRO- 
MACS packagj^ equip ped with a private version of the 
PLUMED plug-irPilSSlj in which the external force 
has been implemented. 

The position zj of the two interfaces is located on-the- 
fly by an algorithm which collects the instantaneous wa- 
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Run 

Face Nu N„ Dp [nm] ^ [nm] a [nm] u/w noi [nm ki 

Aui 

Au2 

Au3 

1.5 

{001} 679 1757 2.5 1.0 0.5 u 2.4 21.0 

3.3 

-^wl 

-^w2 

^w3 

30.25 2.5 

{001} 679 1757 2.5 1.0 0.5 w 28.25 2.5 

26.25 3.5 

NPT 

{001} 679 1757 Ordinary NPT 

Buj 

Bwj 

{110} 625 1533 3.2 1.0 0.2 

w 29.0 1.95 

NPT 

{110} 625 1533 Ordinary NPT 


TABLE I. Simulation settings used for the different MD calculations presented in Sec. |IV[ 


ter density distribution n^{z) in an histogram of bin-size 
5z. Since nv,{z) = 0 inside the crystal volume, we can 
choose a threshold value ni that indicates the transition 
from the liquid phase to the solid. The crystal interfaces 
are then located at n^riz) = ni. Typical values used for 
these parameters are ni = 10nm“^ and Sz = 0.75 A. 

As shown in Fig. to comply with the PBC, is 
applied at both sides of the slab, and the CR consists 
of two layers of liquid at fixed distance from the crys¬ 
tal interfaces. All the C/xMD parameters chosen for the 
presented simulations are based on a preliminary tuning 
performed on a typical system. During this process the 
effective decoupling between the reservoir region and the 
growing crystal is assessed by testing different parame¬ 
ter configurations, and observing the resulting solution 
behavior. All the details of this tuning procedure are 
reported in the Supplementary Material (SM). The rele¬ 
vant calculation settings are listed in Tab. |T] 

IV. CRYSTAL GROWTH AT CONSTANT 
SUPERSATURATION 

We now present the results of the application of the 
C/iMD method to the simulation of urea crystal growth 
in aqueous solution. In our study we have considered 
the growth of the {001} and {110} faces. It is known 
that these two faces are characterized by different growth 
mechanisms: the {001} face undergoes a rough growth 
process, while the {110 } face grows through a birth-and- 

spread mechanisrrP^Iill, 

First, we investigate the {001} face growth process. We 
have performed 2 sets of three simulations, referred to as 
Auj and A^j. In the A^j type simulations restrains 
urea density, while in the A^^j is water to be controlled. 
As reported in Tab. |lj the index j = 1,2,3 runs over 
different target densities. For comparison we have also 
performed an ordinary NPT run of the same system. In 
each simulation we have evaluated the number of solid 
molecules using the Degree Of Crystallinity (DOC) 
variable defined in Ref. [SH 

In Fig. we report the evolution of and of the mole 


fraction x measured in the CR for the A^j and NPT sim¬ 
ulations (an analogous plot for the A^j runs is included 
in the SM) . Let us focus on the behavior of the A^i simu¬ 
lation (upper-left panel). After an initial transient Tt the 
action of F^ stabilizes the mole fraction around x = 0.05. 
As the crystal grows, the reservoir region is gradually de¬ 
pleted and, as a consequence, at Tp = 53 ns the CR mole 
fraction starts drifting to lower values, meaning that the 
C/rMD method is no longer valid. Within the validity 
range (rt < t < rp) increases linearly, as expected 
for a rough growth process at constant supersaturation. 
The behavior during the validity range can be fitted 
with a linear model, the slope of which is the growth 
rate g{ooi} at this particular solution composition. For 
t > Tp the crystal rate decreases, reflecting the change 
in solution chemical potential. An analogous behavior is 
obtained in the A ^2 and Aus simulations, where larger 
supersaturations are enforced, determining a faster crys¬ 
tal growth. This results in a more rapid depletion of 
the reservoir and, as a consequence, in a shorter Tp. In 
contrast with the C/iMD runs, in the NPT simulation 
the solution concentration varies throughout the whole 
dynamics, interfering with crystal growth. In Fig. we 
report the growth rates obtained in the Auj and Awj sim¬ 
ulations, corresponding to different values of x. Remark¬ 
ably, //{ooi} exhibits a linear dependence on x that is char¬ 
acteristic of a rough growth mechanism (see e.g. Ref. I55|l . 
It is rewarding that the results are consistent, irrespective 
from the controlled species. 

In two further sets of simulations we have applied 
C/iMD method to the {110} face growth, either restrain¬ 
ing urea run in Tab. [I]) or water density (Rw)- For 
comparison we have also performed an ordinary NPT 
simulation of the same system. As mentioned before, 
the {110} face growth process exhibits a birth-and-spread 
nature. Because of this mechanism N‘^ evolves in a step¬ 
wise behavior, each step corresponding to the formation 
of a new crystalline layer. The dynamics of iV° for the 
Ru and NPT cases is represented in Fig. In the fig¬ 
ure we have also reported the mole fraction evolution, 
showing that our method is able to maintain a stable so¬ 
lution composition for a significant time. In both Ru and 
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FIG. 4. Auj simulations results compared to the ordinary NPT behavior. The green curves represent the increase of the 
crystal-like molecule number AN‘^ as a function of time. The blue curves represent the solution mole fraction x measured in 
the CR as a function of time. The instantaneous value of x is represented in faded color, while the full-color curves are obtained 
via Exponentially Weighted Moving Average, with characteristic smoothing time of 0.5 ns. The red marks indicate the validity 
time-window of the C/rMD method, namely rt < t < rp. The black dashed lines represent the linear fit of the AN‘^ behavior, 
calculated within the corresponding validity time range as explained in the text. 
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FIG. 5. Growth rates measured in Auj and Awj simulations, 
versus the average GR mole fraction (see Fig. [^. The red 
dashed line is a linear fit of the growth rates. 


Bw runs the CR mole fraction is maintained at approxi¬ 
mately X = 0.062 (all the and results are reported 
in the SM) . 

We now estimate the {110} face growth rate at con¬ 
stant supersaturation 5{iio}- Since the environment so¬ 
lution is not changing, we can assume that after a layer 
growth event the system carries no memory of its past. 


Thus the probability that a new layer is created is in¬ 
dependent from the growth history, and the time inter¬ 
val between two successive events follows an exponential 
distributioiP^. However, in a single run, only few layers 
are created before the solution starts depleting. Thus, in 
order to collect sufficient statistics for the construction 
of the time distribution, we have repeated each B^ and 
B^ simulations hve times. Of course we have considered 
only the events occurring at rt < t < rp. 

We have constructed the cumulative time distribution 
of the observed events, and fitted it with an exponential 
model, as shown in Fig. The statistical significance of 
this ht has been assessed using the Kolmogorov-Smirnov 
testP^, obtaining a p-value of 0.37. From the regres¬ 
sion we have extracted the characteristic occurrence time 
Hiio} = 18.3±2.7 ns, which corresponds to a growth rate 
5{iio} = 1-73 ± 0.25 ns“^. This result can be compared 
to our estimate of the {001} rate ai x = 0.062, that is 
5{ooi} = 6.59 ± 0.97ns“^, obtained from the linear inter¬ 
polation in Fig. 

In Ref. [m we have proposed a model for the growth 
rate of the {hjk} face of urea crystal, in which: 

9{hik} = 90 exp {-l3AG{!,jk}) , 


( 5 ) 
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FIG. 6. Sul simulations results compared to the ordinary 
NPT behavior. The green curves represent the increase of the 
crystal-like molecule number AN'^ as a function of time. The 
blue curves represent the solution mole fraction x measured 
in the CR as a function of time. The instantaneous value 
of X is represented in faded color, while the full-color curves 
are obtained via Exponentially Weighted Moving Average, 
with characteristic smoothing time of 0.5 ns. The red marks 
indicate the validity time-window of the C/rMD method. 



FIG. 7. Gumulative time distribution of the layer growth 
events a.t x = 0.062 {Bu and Bw simulations). The red curve 
is the fitted Cumulative Distribution Function, equal to y{t) = 
1—exp(—t/r{iio}) for exponentially distributed time intervals. 

where go is a characteristic diffusion limited growth rate, 
^G^hjk} is the free-energy barrier associated to the cre¬ 
ation of an extra {hjk} layer and = ^bT. The model 
has been used in Ref.lT^to predict the urea crystal habits, 
by evaluating the velocity ratio between different crystal 
faces: 

9 {hj k} / 9 {Imn} tiXp( (6) 

We underline that this ratio depends only on the free- 


energy barriers, which are static properties of the sys¬ 
tem. Using the AG{ooi} and AGmoi evaluated via Well- 
Tempered (WT) metadynamic^^in Ref.Hni Eq. (|^ gives 
5{ooi}/5{iio} = 3-61 at a; = 0.062. This result is consis¬ 
tent with the ratio obtained from our dynamical esti¬ 
mates, that is (/{ooil/llliio} = 3.8 ± 0.8. 

The birth-and-spread nature of the {110} face growth 
mechanism is determined by non-negligible free-energy 
barriers AGi associated to the incorporation of new so¬ 
lute molecules in the growing crystal. In the following 
we shall derive such free-energy barriers from the C/rMD 
trajectories obtained from the and simulations. 

The growth dynamics, see e.g. Fig. shows that the 
system evolves through a sequence of metastable states 
i = 1.. .n, each characterised by an average number of 
molecules in the crystal state Therefore, the prob¬ 

ability distribution P(A°), computed averaging over the 
ensemble of i?u and B^ trajectories, will exhibit an alter¬ 
nation of minima and maxima, the former representing 
the metastable states at and the latter 

representing the faster layer growth transitions. This is 
in agreement with Fig. which represents the behavior 
of: 

G(iV°) =-^lnP(lW). (7) 

From the function G{N‘^) we can evaluate the free-energy 
barriers associated to the layer growth transitions. Since 
the lifetime of each metastable state is such that the 
sampling in the vicinity of a minimum can be consid¬ 
ered ergodic, the relative free energy associated to the 
fluctuations within the i-th basin can be written as: 

= -i In = G(JV=) - G(JVG„,,) 

( 8 ) 

According to Eq. (|^, the free-energy barrier associated to 
the creation of a layer from a metastable state is AGi = 
G(Afmax,J - as indicated in Fig.[^ 

The barriers resulting from our C/xMD simulations ex¬ 
hibit very similar heights, showing that successive growth 
events obtained in the molecular model are equivalent. 
This provides further confirmation that the action of 
C/iMD maintains the growing crystal environment at 
a constant chemical potential. As shown in Fig. the 
data can be fitted using a sinusoid, obtaining a remark¬ 
able agreement. The free-energy barrier associated to 
the sampled growth events can be extracted from the 
fitted curve, obtaining AG = 2.63 k^T, which is in sub¬ 
stantial agreement with the corresponding barrier AG = 
2.0±1.0 k^T, computed in the model proposed in Ref.lHl 


V. CONCLUSIONS AND PERSPECTIVES 

In this work we have presented the C/rMD method, 
which addresses the concentration depletion problems in 
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FIG. 8. Behavior of G{N‘^) (see Eq. Q) evaluated using the 
probability distribution P{N^), extracted from the and 
Bw simulations sampling. The red dashed line represents a 
sinusoidal fit of the numerical result. 


the MD study of finite systems. C//MD applies an exter¬ 
nal force that controls the solution composition within 
a region of interest of the system, while the remaining 
volume is used as a molecular reservoir. 

C//MD has been implemented and tested for the rele¬ 
vant case of urea crystal growth in aqueous solution. For 
the considered systems the method is capable of main¬ 
taining the solution environment of the growing crystal at 
constant composition, up to times of the order of 10-f-100 
ns. This has allowed us to estimate the {001} and {110} 
faces growth rates in a constant supersaturation envi¬ 
ronment. Such kind of results are not accessible with 
ordinary NPT dynamics, due to the intrinsic coupling 
between the crystal size and the solution composition. 
The retrieved results are also in remarkable agreement 
with the predictions of Ref. [T^ which are based on the 
combination of a theoretical free energy model and WT 
metadynamics calculations. The evaluation of the free- 
energy barriers associated to the growth of consecutive 
{110} layers has shown that the events are thermody¬ 
namically equivalent, providing further support to the 
hypothesis that C^MD can establish a stationary growth 
regime. 

The proposed scheme can be applied to a variety of MD 
problems, such as crystal nucleation, electro-chemistry, 
surfactants adsorption. However the external force has 
to be properly re-defined to comply with the character¬ 
istic features of these systems. The timescale of validity 
of C/xMD is related to the fixed number of molecules in¬ 
volved in the simulation. A possible future development 
can be th e com bination of C/rMD method with adaptive 
resolutiorP^^ or particle insertion techniqueJ^, which 
would result in a considerable extension of the validity 
range of the method. 
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In the following we report further details about the results presented in the Main Submission (MS) text. This 
supplementary material is organized in two sections: Sec. [^presents a preliminary study that has been performed to 
tune the C/rMD method for the typical system considered in our calculations. In Sec. |SII| we gather supplementary 
information concerning the simulation results presented in the MS. 


SI. PRELIMINARY STUDY 

In this section we report the results of a preliminary 
study on the application of C/rMD method to crystal 
growth. This study aims at finding convenient parameter 
settings for the typical urea-water systems simulated in 
our work. 

It has to be noted that C/rMD reliability is related to 
system properties such as size, initial solution concen¬ 
tration and liquid diffusivity. For this reason, a tuning 
process becomes necessary each time the method is ap¬ 
plied to a new system. 

Let us first summarize the C/xMD parameters: 

• Dp, the distance of the force center Zp from the 
solid-liquid interface zj. 

• the CR position and size, defined by ^ and tr. The 
former is the distance between the inner CR bound¬ 
ary and zi, while the latter is the distance between 
the outer CR boundary and Zp (see Fig. 2 in the 
MS). 

• the force constant ki (see MS Eq. (1)), 

• the characteristic length lu of MS Eq. (4), 

• the target concentration ngi, either for urea or wa¬ 
ter species. 

The system used for this study is the {001} growth con¬ 
figuration, all the relevant parameters are summarized in 
Tab. In particular, the following analysis is focused 
on tuning of Dp, a and while constant values have 
been assigned to riuO; C w. In order to keep a limited 
growth rate, and thus a long C/iMD validity time, n^o 
has been set to 1.5 nm“^. The choice for the ^ and w 
values will be discussed later on. 

In a first set of 4 MD runs, which we will refer to as Xj, 
we have tested the effect of varying Dp and a together, 
so that the force is applied at different distances from 
the crystal interface while the CR maintains the same 
size and position (referred to the crystal). 

In Eig. we report the evolution of the CR urea den¬ 
sity It can be noted that, in each simulation, the 

C^MD maintains a stable close to the target value 
riuo, for at least 52 ns {X 3 and X 4 ). In the Xi and X 2 
cases, for which is applied closer to the CR, the va¬ 
lidity time is longer. This is because a larger volume is 


assigned to the reservoir region, and the control of the 
solution concentration is more effective. 


To look at the quantitative differences between Eig.|Sl| 
curves, in Tab. |SII| we report the average CR density 
and standard deviation for each run, evaluated over 52 
ns. During this lapse of time all the simulations enforce 
an average density close to the target value nuo, and a 
correlation between the density fluctuations and a (or, 
alternatively. Dp) emerges. That is because a larger dis¬ 
tance between the force center and the CR results in a 
slower response of F^ to the density variations, and thus 
in a weaker control of 


In the Xj runs we have also monitored the mole frac¬ 
tion at different distances from the crystal. In Eig.|S2|we 
report the resulting mole fraction profile for each simula¬ 
tion, averaged in time. The qualitative behavior of this 
profile is consistent with the trend described in Sec. 2 of 
the MS and, in particular, we underline that the choice 
of ^ = 1.0 nm is appropriate to identify the TR of our 
system. 


Fig. [S^ shows that, for the X 2 , X 3 and X 4 cases, the 
application of C/rMD method determines a flat profile 
within the CR, and no correlation with Dp is present. 
Conversely, in the Xi case the vicinity of the force cen¬ 
ter to the CR has a direct contribution on its concen¬ 
tration. This should be avoided, since it prevents the 
decoupling of the reservoir region from the environment 
of the crystal. To support this statement we compare the 
time evolution of the x profile for the Xi and X 2 cases. 
To this purpose we have evaluated the time average of 
the X profiles over 4 consecutive windows of time, so that 
the slow composition dynamics is visible. Let us focus on 
the behavior of x within the CR region: while in the Xi 
case (represented in Fig. S3i) the profile varies through¬ 
out the whole simulation, in the X 2 case (Fig. S3 o) the 
mole fraction remains flat for the first three windows of 
time, up to the limit of validity of C/rMD. 


These results suggest that a layer of thickness a be¬ 
tween the force center and the CR is necessary to decou¬ 
ple the reservoir from the crystal environment. However, 
the choice of tr also affects the promptness of in bal¬ 
ancing the CR composition variations. Given this, the 
configuration used in X 2 case seems to be the best com¬ 
promise. We underline that the choice of a is also related 
to the characteristic length w, which indicates the local¬ 
ization of around the force center. In our calculations 
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Run 


Nu Ny, Df [nm] ^ [nm] a [nm] w [nm] nou [nm fcu 


Xi 


2.0 

0.0 


X 2 

679 1757 

2.5 

1.0 0.5 0.2 1.5 

21.0 

A 3 


3.0 

1.0 


A 4 


3.5 

1.5 


NPT 2 V 

698 5125 


Ordinary NPT, double size 


Fi 




35.0 

F 2 

679 1757 

2.5 

1.0 0.5 0.2 1.5 

21.0 

F 3 




7.0 

Yi 




3.5 


TABLE SI. Simulation settings used for the MD simulations reported in this section. 



FIG. SI. Trajectories of the CR urea density for the Xj simulations. The instantaneous value is represented in faded color, 
while the full-color curves are obtained via Exponentially Weighted Moving Average, with characteristic smoothing time of 0.5 
ns. Each curve is referred to the target density nou = 1.5nm“® (dot-dashed axes). The dashed vertical line indicates t — 52 
ns, the time at which solution depletion becomes significant for X3 and X4. 


Run 

(njf^) - nuo 

^71 

Ai 

-2.0 X 10"® 

0.24 

A 2 

3.1 X 10“^ 

0.26 

As 

3.8 X 10~2 

0.30 

A 4 

4.9 X 10“^ 

0.33 


TABLE SII. Average CR urea density (referred to Uuo = 1.5) 
and standard deviation Sn = \/((^u^)^) “ of fh® 

simulations. The density, in units of nm“®, is averaged over 52 
ns of dynamics, in order to consider only the validity regime 
of the method. 


we have chosen w = 0.2-L 0.3 which, according to our ex¬ 
perience, determines an efficient, but still local, action on 
the molecules. 

In a further set of 4 simulations, indicated in Tab. 
as Yj, we have performed C/rMD with different values of 
the force constant starting from the X 2 settings. As a 


term of comparison we have also performed an ordinary 
NPT simulation, named NPT 2 y, in which the crystal 
slab is immersed in a larger volume of liquid, almost twice 
the original size, with a solution composition close to that 
enforced by in the Yj runs. 

In Fig. we report the evolution of the CR urea con¬ 
centration n™ in the Yj and NPT 2 y simulations. First 
we note that the NPT 2 V system does not experience a 
significant solution depletion up to 20 ns. Thus, within 
this time range, we can consider the NPT 2 y run as a 
good approximation to a macroscopic system. 


If we look to the Yj trajectories we see that the 
dynamics is affected by the choice of k. In the Yi case 
(fcu = 35.0nm^ kJ/mol) fluctuates steadily around 
riuo, while for smaller fcu the density is less stable. This is 
quantitatively confirmed by Tab. |SIII where we report 
the corresponding riu^ mean and standard deviations, 
up to 20 ns. As fc is reduced, the standard deviation 
increases, while the average density distances the target 
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FIG. S2. Local mole fraction profile as a function of the 
distance from the solid-liquid interface. The results of the 
Xj simulations, averaged over r = 41 ns, are displayed. Each 
local value is evaluated in 0.25 nm thick layers, symmetrically 
located at both sides of the crystal. The vertical colored lines 
highlight the force center position except for the X 4 , case, in 
which Df falls outside the plot range. The error bars are 
evaluated with the re-blocking techniqu^^. Because of the 
long correlation times the error estimate is not completely 
reliable for D >2 nm. 



0.5 1.0 1.5 2.0 2.5 3.0 


D [nm] 


Run 

(^u ) ^uO 

Fi 

1.5 X 10"^ 0.24 

F 2 

5.5 X 10"^ 0.26 

F 3 

7.5 X 10"^ 0.29 

Yi 

1.8 X 10"^ 0.34 

NPT 2 V 

1.6 X 10"^ 0.34 


TABLE Sill. Average CR urea density (referred to Uuo = 1.5) 
and standard deviation of the Yj and NPT 2 V simulations. 
The density, in units of nm“®, is averaged over 20 ns of dy¬ 
namics. 


sity fluctuations in the realistic NPT dynamics are more 
similar to the results of I 4 simulation, the one with the 
smallest k. This raises the question whether this differ¬ 
ence in composition sampling affects the crystal growth. 

To answer this question in Fig. |S5| we report the lo¬ 
cal standard deviation s„(ZI) of . While at distances 
D > 1.5 nm from the interface the Yj runs exhibit a 
different s„ behavior compared to the NPT 2 y case, this 
difference is much mitigated as we get closer to the solid- 
liquid interface. This analysis suggests that, for this sys¬ 
tem, the choice of does not affect significantly the 
composition fluctuations of the solution in contact with 
the crystal. As a result it is preferable to choose a large 
fcu, which entails a more effective control on the CR. 

As a last remark we underline that in this study 
was applied exclusively to the urea specie. However the 
results are still useful when the restraining force acts on 
the water population, provided that k is rescaled accord¬ 
ing to the typical riw fluctuations. 


Sll. CRYSTAL GROWTH CALCULATIONS 

In this section we report the additional material con¬ 
cerning the MD calculations presented in Sec. 4 of the 
MS. We organize the information in two sub-sections, 
concerning the growth processes of the { 001 } and { 110 } 
faces of urea. 


FIG. S3. Temporal evolution of the average mole fraction 
profile for the Xi (a) and X 2 (b) cases. Each curve represents 
the average of x{D) over a different time window, as indicated 
in the legend. The black arrows highlight the temporal evo- 
Intion of the composition in the reservoir. The error bars are 
evaluated with the re-blocking techniqu^^. Because of the 
long correlation times the error estimate is not completely 
reliable for D >2 nm. 


value riuo- That is because a larger fcu determines a more 
intense response to the density fluctuations, and thus 
a more effective control on the flux of urea molecules from 
the reservoir to the CR. 

From Tab. |SIII| we also note that, while a larger k 
allows a better control of the CR composition, the den- 


A. {001} face growth. 


For the study of {001} face growth we have performed 
two sets of simulations, named A^j and (see Tab. I 
of MS). In the first one the urea population is controlled, 
while in the second water is controlled. The results of the 
Auj set are reported in the MS, here we show the results 


of the A^j runs, in Fig. S 6 compared to the ordinary 


NPT dynamics. Also in this set of simulations the C/rMD 
method succeeds in establishing a linear growth regime, 
even though the control of x is less rigid and allows larger 
fluctuations. 

The dependence of the growth rates extracted from 
both the Auj and A^j simulations on the corresponding 
solution mole fraction has been assessed with a linear fit. 
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FIG. S4. Trajectories of the CR urea density for the Yj and NPT 2 V simulations. The instantaneous value is represented in 
faded color, while the full-color curves are obtained via Exponentially Weighted Moving Average, with characteristic smoothing 
time of 0.5 ns. Each curve is referred to nou = 1.5nm“^ (dot-dashed axes). The dashed vertical line indicates t = 20 ns, when 
solution depletion becomes significant in the NPT 2 V run. 



FIG. S5. Logarithmic plot of the standard deviation, for 
the Yj and NPT 2 V runs. Each value is measured by detecting 
nu within a 0.25 nm thick layer placed at a distance D from 
the solid-liquid interface. The dashed vertical lines indicate 
the GR boundaries. 


represented in Fig. 5 of the MS. The parameters resulting 
from the interpolation are the following: 

g =-0.83 ± 0.75 ns"\ 
m = 119.8 ±9.6 ns“^, 

where the linear model \s y = mx ± q. 


5 repetitions of the same simulations, in order to collect 
a relevant statistics for the presented analysis. Moreover, 
the C/rMD settings have been designed so that both B^j 
and i?wj simulations reach analogous growth conditions. 
For this reason, in our analysis, we have considered the 
results of the and B^j runs as equivalent. In the 
multiple plot of Fig. [^we report the resulting dynamics 
of each simulation. 

In the MS we have reported the estimate of the {110} 
face growth rate, obtained from the trajectories displayed 
in Fig. |S7| As explained in the text, under constant su¬ 
persaturation conditions we can assume that the prob¬ 
ability of growth events occurrence follows a Poisson 
distribution. Therefore the waiting time between two 
consecutive events is exponentially distributed, and the 
growth rate can be extracted from this distribution. To 
this purpose the time intervals between two consecutive 
growth events have been collected. We have considered 
only those events occurring within the validity time of 
C/rMD, during which the growth could be considered at 
constant supersaturation. We have also discarded the 
waiting time preceding the first event, since it is affected 
by the initial transient time, during which the solution is 
brought at the desired composition. 

From the resulting statistics we have then constructed 
the cumulative time distribution, and fitted it with 
the cumulative distribution associated to an exponential 
probability density: 


B. {110} face growth. 


y{t) = 1 - exp(-t/T{iio}). (SI) 


For the study of {110} face growth the B^j and B^j 
simulation sets have been performed, the relative param¬ 
eters are shown in Tab. I of MS. Each set is constituted by 


From the fit we have extracted ''"{no} = 18.3 ± 2.7 ns 
(as shown in Fig. 7 of the MS). The corresponding er¬ 
ror has been obtained with the bootstrap techniques!! 
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t [ns] t [ns] 

FIG. S6. Av,j simulations results compared to the ordinary NPT behavior. The green curves represent the increase of the 
crystal-like molecule number AN‘^ as a function of time. The blue curves represent the solution mole fraction x measured in 
the CR as a function of time. The instantaneous value of x is represented in faded color, while the full-color curves are obtained 
via Exponentially Weighted Moving Average, with characteristic smoothing time of 0.5 ns. The red marks indicate the validity 
time-window of the C/rMD method, namely rt < t < rp. The black dashed lines represent the linear fit of the AN‘^ behavior, 
calculated within the corresponding validity time range. 


The growth rate has been estimated by calculating the 
average increase associated to the growth events 
(31.6±0.3), and dividing it by T{iio}, to obtain the result 
reported in the MS. 


no. 1, 1989. 

[S2]B. Efron, “Bootstrap methods: Another look at the jackknife,' 
The Annals of Statistics, vol. 7, no. 1, pp. pp. 1-26, 1979. 
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FIG. S7. Bnj and B^j simulations results. The green curves represent the increase of the crystal-like molecule number 
as a function of time. The blue curves represent the solution mole fraction x measured in the CR as a function of time. The 
instantaneous value of x is represented in faded color, while the full-color curves are obtained via Exponentially Weighted 
Moving Average, with characteristic smoothing time of 0.5 ns. The horizontal dashed line shows the average of x, evaluated 
over the validity time-window of the C/rMD method, indicated by red marks. 









































































































